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Abstract 



Predictions of the future state of the Earth's atmosphere suffer from the consequences of chaos: numerical weather forecast models 
quickly diverge from observations as uncertainty in the initial state is amplified by nonlinearity. One measure of the utility of a 
forecast is its shadowing time, informally given by the period of time for which the forecast is a reasonable description of reality. 
The present work uses the Lorenz '96 coupled system, a simplified nonlinear model of atmospheric dynamics, to extend a recently 
developed technique for lengthening the shadowing time of a dynamical system. Ensemble forecasting is used to make forecasts 
with and without inflation, a method whereby the ensemble is regularly expanded artificially along dimensions whose uncertainty 
is contracting. The first goal of this work is to compare model forecasts, with and without inflation, to a true trajectory created 
by integrating a modified version of the same model. The second goal is to establish whether inflation can increase the maximum 
,— ' shadowing time for a single optimal member of the ensemble. In the second experiment the true trajectory is known a priori, and 

| | only the closest ensemble members are retained at each time step, a technique known as stalking. Finally, a targeted inflation is 

introduced to both techniques to reduce the number of instances in which inflation occurs in directions likely to be incommensurate 
Q-iwith the true trajectory. Results varied for inflation, with success dependent upon the experimental design parameters (e.g. size 
(— i of state space, inflation amount). However, a more targeted inflation successfully reduced the number of forecast degradations 
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without significantly reducing the number of forecast improvements. Utilized appropriately, inflation has the potential to improve 
predictions of the future state of atmospheric phenomena, as well as other physical systems. 

Keywords: shadowing, stalking, ensemble forecasting, inflation 



1. Introduction 

Society is often dependent upon the ability of scientists to 
accurately forecast the future state of chaotic physical systems, 
such as the weather. In extreme cases, atmospheric scientists 
are asked to anticipate natural disasters with ample time for the 
public to adequately prepare. To predict the likelihood of at- 
mospheric phenomena, meteorologists estimate Earth's future 
atmospheric state using sophisticated computer simulations. In 
general, these mathematical models struggle to predict the spe- 
cific behavior of chaotic physical systems due to (1) uncertainty 
in the initial state, (2) chaos, and (3) model error yj]. 

The initial state is an estimate of the condition of the at- 
mosphere at the beginning of a forecast. Although weather 
monitoring devices observe much of the planet, for many ar- 
eas (e.g. southern hemisphere oceans, the upper atmosphere) 
the measurements are spatially sparse and typically three orders 
of magnitude less in number than the degrees of freedom in a 
global weather model. Therefore, meteorologists must combine 
observations with prior forecasts to estimate initial conditions 
in a process called data assimilation. 
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Since perfect knowledge of atmospheric tendencies is un- 
achievable, a great deal of recent research has focused on data 
assimilation, the process by which observations are combined 
with model predictions to give the best possible initial state, 
or analysis |2j-|6J]. It is the analysis that is typically used as a 
proxy for the true state of the atmosphere at any time in the 
past. Although the analysis has inherent uncertainty from lack 
of perfect observations and an imperfect model, the modeler's 
goal is to create a forecast from this given state that remains 
reasonable (i.e. close to observations) for the longest duration 
of time. In the context of weather predictions, forecasters strive 
to be accurate for two weeks, the limit imposed by chaos. To- 
day, five-day forecasts are as good as three-day forecasts were 
30 years ago 01. 

Despite these efforts, the mischaracterization of Earth's at- 
mospheric state introduces uncertainty. In addition, the atmo- 
sphere is a chaotic system, causing these small uncertainties in 
the initial state to grow exponentially in time. Finally, since the 
model itself is not a perfect representation of reality, inaccurate 
changes forecast by the model, namely model error, compound 
the uncertainty and lead forecasts to diverge quickly from ob- 
servations. As Lorenz noted in 1965, the limit of predictability 
of the atmosphere is about two weeks, even with nearly perfect 
knowledge of the current state [8]. 

To account for initial state uncertainty, an accepted technique 
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is to use ensemble forecasting, where a collection of perturba- 
tions to the best guess are chosen randomly from a distribution 
that reflects the measurement uncertainty and system dynamics 
near the given initial state [9]. Each ensemble member is fore- 
cast forward in time, yielding a collection of final states. This 
collection roughly represents the probability distribution of the 
model's forecast with associated uncertainty. For example, if 
60% of the ensemble members predict rain, the forecaster as- 
signs a 60% chance of rain. Given the limitations, the modeler's 
goal can be stated as trying to keep some ensemble members 
close to the observed truth for as long as possible, analogous to 
finding a shadowing trajectory for Earth's atmosphere. 

Traditional mathematical shadowing theory was developed 
for hyperbolic systems, based on the Shadowing Lemma of 
Anosov II 1 OTI and formalized by Bowen 111 ill . Given a pseudo- 
trajectory of a model (i.e. one very close to an actual one-step 
trajectory), this lemma establishes the existence of a true tra- 
jectory that remains close for an arbitrary period of time. Later 
research has extended the lemma for a wide variety of hyper- 
bolic systems (e.g. Ill2l4l4ll ). For these systems, the number 
of expanding directions remains constant, i.e. the number of 
Lyapunov exponents greater than zero is constant throughout 
the state-space, and small perturbations in stable directions de- 
cay exponentially in time. However, most physical systems 
(e.g. Earth's atmosphere) are non-hyperbolic, i.e. the num- 
ber of positive Lyapunov exponents fluctuates along a trajec- 
tory through state space. For these systems, there does not exist 
a trajectory that shadows the truth for an arbitrarily long time 
1 15141711 . More recent work has been focused on finding shad- 
owing trajectories for low-dimensional non-hyperbolic systems 
(e.g. general ordinary differential equations 118lll9ll . the driven 
pendulum [20], the Henon map 12111 '). 

For the purposes of the present research, the time period for 
which a particular forecast is an accurate representation of real- 
ity is referred to as the shadowing time. Danforth and Yorke 
1 22] proposed an aggressive approach to shadowing called 
stalking to increase the shadowing time for a given forecast. For 
a system with 71-degrees of freedom, an ^-dimensional sphere 
was used to encompass the initial ensemble members. Through- 
out the length of the forecast, these members experience expan- 
sion away from each other in some dimensions and contraction 
towards each other in others, depending on the local finite-time 
Lyapunov exponent in each dimension. The ensemble forecast 
can be approximated by an n-dimensional ellipsoid for a short 
time. 

The idea of stalking is to artificially impose some additional 
uncertainty in the contracting directions (i.e. those with a neg- 
ative local finite-time Lyapunov exponent) at periodic inter- 
vals throughout the forecast. This is accomplished by inflating 
the ellipsoid along axes parallel to the contracting dimensions. 
Thus, if the true state of the system happens to suddenly expand 
along a previously contracting direction, as happens in systems 
exhibiting unstable dimension variability lE3ll . some ensemble 
members will remain relatively close to the true state. Stalking 
is not currently used in weather forecasting, but an alternative 
version called variance inflation is used in ensemble-based data 
assimilation to ensure the state estimation algorithm does not 



put too much faith in the model forecasts and ignore observa- 
tions when creating the analysis 112.411 . 

For the purpose of this paper, the terms "shadowing" and 
"stalking" are reserved for experiments where the truth is 
known a priori, and the goal is to identify an optimal ensem- 
ble member (section V). Using the technique described above 
to create forecasts, without prior knowledge of the truth, is here 
denoted forecasting with and without inflation (section IV). In 
section II, we present the non-hyperbolic '96 Lorenz system 
1 25] and the model, and examine the effect of changing basic 
parameters. In section III, the process of creating the initial en- 
semble and making forecasts is explained. In section IV, the 
inflation technique is introduced, and forecasting results (with 
and without inflation) are presented. In section V, we present 
the results from stalking experiments. In section VI, we use tar- 
geted inflation to isolate the instances in which inflation is most 
likely to improve a forecast. Finally, the findings are discussed 
with respect to operational weather prediction in section VII. 

2. The Lorenz 1996 System 

In numerical weather prediction (NWP), the true evolution 
of the Earth's atmosphere will ideally be a plausible member of 
an ensemble forecast [9]. In our context, ensemble forecasting 
is used to model the trajectory z fl , the "truth," of some meteo- 
rological quantity (e.g. temperature, pressure). In the present 
study, both the truth and the forecast were created using ver- 
sions of a simplified nonlinear model given by Lorenz (1996) 
to represent the atmospheric behavior at / equally spaced loca- 
tions on a given latitude circle 02511 . This system has been used 
to represent weather related dynamics in many previous studies 
(e.g. I26T-280. The governing first-order differential equations 
are given by 12511 : 
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for i = 1,2, ...,/, j = 1,2, ...,//, andra = (7+1)7. Although n 
represents the full state space size, we are primarily interested 
in the values of the slow variables x,. Thus, for the purposes of 
this study we consider I to be the dimensionality of the system. 

The values jc, represent slowly changing meteorological 
quantities whose dynamics are described by Eq. ([Tal l. Since the 
set of xi corresponds to locations along a single latitude circle, 
the subscripts i and j are defined to be in a cyclic chain. That is, 
we define X-\ = Xr-i, xq - Xj, and x\ - Xj+\, and similarly for 
j. Each xi is then coupled to 7 quickly changing, small ampli- 
tude variables whose behavior is governed by Eq. (TTbV For our 
experiments, we set I — 4, 5, 6 and 7=16 yielding n = 68, 85, 



2 




Figure 1: Model Schematic t26ll . The 7=8 slow variables Xj can be thought of 
as locations along a latitude circle. Each variable has J = 4 corresponding yt 
fast variables for a total of n = 40 degrees of freedom. 



and 102 variables. A schematic is shown in Fig.[T]where eight 
slow variables (x,) are each coupled to four fast variables (yj), 
reflecting a total of n — 40 degrees of freedom. Note that the 
dynamics of each x, are dictated by neighboring x variables and 
the corresponding set of coupled yj variables. 

As a low-dimensional model of atmospheric dynamics, this 
system is ideal for several reasons; primary among them being 
the ability to tune the relative levels of nonlinearity, coupling of 
timescales, and spatial degrees of freedom. The nonlinear terms 
in Eq. dTab are meant to represent advection, and conserve the 
total energy of the system. The linear term signifies a loss of 
energy either through mechanical or thermal dissipation. Ex- 
ternal forcing F is then added to prevent the total energy from 
decaying completely. For all experiments we set F = 14. Con- 
sistent with the literature, we set c = 10 and b — 10, which 
forces the fast variables to oscillate 10 times quicker than the 
slow variables Il7ll26l 12811 . The coupling parameter is set to a 
value of h — 1 for the system truth, and h = 0.5 for the model. 
Note that one time unit corresponds to five days, the dissipative 
decay time I29I1 . 

Integration of the differential equations is completed using 
the fourth order Runge-Kutta method with a time step of 0.01. 
Rigorous shadowing attempts could be made using far more 
advanced methods of integration, with much smaller time steps 
1 30]. However, for the purpose of this study of short forecasts, 
the difference is negligible. 

2.1. Dynamics 

The system dynamics can be observed by integrating 
Eqs. ( [Tal l and ( fTbb with h = 1, hereafter referred to as the sys- 
tem. In Fig. |2l a time series for a particular longitudinal pro- 
file (/ = 40, J = 16) is shown after variable X\j, is initially 




(b) Days 50 - 55 

Figure 2: A time series of the system shows the effect of a five-unit perturbation. 
Panel (a) shows days 0-5 and panel (b) shows days 50 - 55. Perturbations 
appear to move primarily to the east (right). After five days, perturbation effects 
can be observed at over half of all the locations (/ = 40, J = 16). 



perturbed by five units. Profiles are recorded at 12-h intervals 
for days 0-5 (panel a) and days 50 - 55 (panel b). Advec- 
tion is apparent as the energy from the perturbation is observed 
halfway around the latitudinal circle by day five. Lorenz and 
Emanuel [29] calculated a perturbation growth rate (doubling 
time) of approximately two days for this model, which agrees 
with trends in larger atmospheric models where errors double 
in roughly two days ll3lll . However, growth rates over limited 
time intervals as in Fig. |2(b)| can differ greatly depending on / 
and J. 

Adjusting the parameters of the system, and examining the 
time series at a single latitudinal site, the interdependence of 
the slow and fast variable coupling is observed. For example, 
with 1-6 and J — 16, a "regular" oscillatory pattern occurs as 
an energy equilibrium is achieved between external forcing and 
dissipation. In this work, we will use the term "regular" to refer 
to these quasi-periodic patterns. However, as the significance 
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of the fast variables is reduced by lowering the coupling param- 
eter h, the time series becomes more complex. Less energy is 
able to dissipate from slow variables to fast, leading to greater 
irregularity in x,-. 

Similarly, the regularity of a time series can be adjusted by 
varying the number of slow or fast variables. As illustrated in 
Fig. [3] (across columns), increasing J creates a more regular 
system for all values of 7. For 7 = 6, the system is fairly regular 
for both J = 16 and J = 40. However, with only a few fast 
variables active (e.g. J = 8), oscillations become increasingly 
aperiodic. For 7 = 8, additional fast variables are required to 
achieve the consistent pattern; they allow the system's energy 
to be more evenly distributed, and nonlinearities in the fast vari- 
ables have less impact on the stability of the slow variables. 

Adjusting the number of slow variables greatly alters the 
qualitative dynamics of the time series for each particular x, 
(Fig.[3j across rows). For 7 = 8, only 7 = 4 exhibits regularity. 
However, with additional fast variables (J = 16), the 7 = 6 sys- 
tem exhibits the greatest regularity. For 7 = 4, the slow variable 
relative maxima are more pronounced than the relative minima, 
and the system appears to be quasi-periodic with a dominant 
frequency of roughly ten days. For 7 = 6, the symmetrical pat- 
tern emerges with a frequency of roughly five days. This pattern 
is similar when J is increased to 40, but with decreased ampli- 
tude. However, for 7 = 4 and 7 = 8 with J = 40, the pattern 
being repeated is more complex. 

We make the preceding observations to motivate the choice 
of Lorenz '96 as an appropriate testbed with which to perform 
ensemble forecasting experiments. The dynamics are quite sen- 
sitive to the system parameters, and we will describe forecast 
results, which vary depending on the underlying dynamics. 

2.2. Coupling 

Each Xi in the system is coupled to J small-amplitude fast 
variables yj, and the strength of this coupling is controlled by 
the parameter h. As mentioned above, the system is defined by 
setting h = 1 . However, the model is rendered imperfect by 
dampening the effect of the fast variables in Eq. ( TTbl by setting 
h = 0.5, hereafter referred to as the model. A 20 day time se- 
ries for a single x, with its corresponding y/s is shown in Fig. [4] 
(7 = 8, J = 4). In the top frame we compare the time series 
of x\ for integrations of the same initial condition by both the 
system and the model. By day 20 the x\ values differ by the 
climatological span of x\ . In addition to increasing the ampli- 
tude of slow variable oscillations, reducing the coupling h has 
the effect of creating a less regular time series. Note that the 
fast variables (shown for the system) vary with amplitudes on 
the order of 10% of the slow variables. 

3. Methods 

From a given initial condition, the trajectory of the truth (z") 
is created by integrating the system in Eqs.tfTal and ([TBI . Two- 
and three-dimensional views of orbit segments of this attractor 
with 7 = 4, J = 16 are shown in Fig. |5(a)| and |5(b)| As men- 
tioned above, the forecast for each ensemble member is then 




Time (days) 

Figure 4: A 20 day time series for x\ and 3*1,2,3,4 using / = 8,/ = 4. For 
the top frame, the solid line represents the time series for X\ with fast variable 
coupling set at h = 1 (the system). The dashed line has fast variable coupling 
reduced to h = 0.5 (the model). The bottom four frames are the fast variable 
time series for the system yi.2,3.4 coupled to x\ , as illustrated in Fig.[T] 

created by reducing the weight of the fast variables by 50% in 
the model, relative to the system. Two- and three-dimensional 
slices for the model can be seen in Fig. |5(c)| and |5(d)| This 
particular experimental design was chosen as it is typical for 
global atmospheric models to attempt to parameterize sub-grid 
scale behavior (e.g. for phenomena occurring on a finer tempo- 
ral/spatial scale). 

3.1. Ensemble creation 

The experimental design is inspired by and follows that of 
Wilks B26I1 . First, long integrations of the system on random- 
ized initial values were performed to establish the shape of the 
attractor. A set of 500 different 7-dimensional points were then 
chosen at 250-day intervals. This spacing was chosen to ensure 
that the initial states sampled different regions of the system 
attractor, and neighboring states were uncorrected. A 'true' 
trajectory z" from each of these states was determined using a 
50-day integration of the system. The goal is to use the model 
to shadow each of these individual true trajectories with an en- 
semble of 20 members. Care was taken to establish an ensemble 
that was a dynamically realistic sample of state space near the 
initial state as follows. 

At each of the 500 initial states, an 7-dimensional hyper- 
sphere was constructed encompassing 100 neighboring states 
from the system attractor, identified during a very long system 
integration. A neighboring state is defined to be one within 5% 
of the span of x, in the i th dimension of the attractor. For each 
hypersphere, the covariance of the 100 neighboring states is cal- 
culated, yielding a 7 x 7 matrix C. The matrix values are then 
scaled to ensure that the average ensemble standard deviation is 
5% of the climatological span of the system attractor. Thus, for 
each of the 500 hyperspheres, 20 initial ensemble members are 
drawn from the distribution: 
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Figure 3: A 50 day time series for x$ illustrates the dependence of the system stability on the number of slow (7) and fast (7) variables. The number of slow variables 
varies from 7 = 4 (a-c) to 7 = 6 (d-f) to 7 = 8 (g-i). The fast variables vary from 7 = 8 (left column) to 7 = 16 (middle column) to 7 = 40 (right column). Generally, 
increasing 7 leads to increased regularity, while changes in 7 have an impact that is dependent upon the magnitude of 7. 
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(a) Two-dimensional view of orbit segment of system attractor (b) Three-dimensional view of orbit segment of system attractor 




(c) Two-dimensional view of model attractor (d) Three-dimensional view of model attractor 



Figure 5: System (a-b) and model (c-d) attractor. Panels (a) and (c) show two-dimensional views looking at xi vs X3. Panels (b) and (d) show three-dimensional 
views using xi,X2, and X4 (I = 4, / = 16). 
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where A is the average eigenvalue of C and t is the standard 
deviation of x, observed in the attractor. 

First, a control state is picked within each hypersphere by 
adding appropriately distributed random noise to the / slow 
variables of the truth as follows: 



z f Ji, 1) = z a M, 1) + 



lt y(i, i), 



(3) 



for i = 1,2,...,/. The vector y is /-dimensional, consisting 
of random entries from a Gaussian distribution. The Cholesky 
decomposition is used to calculate the square root of C lnlt [32]. 
Note that if C; n it was not symmetric positive definite (due to 
the finite sample in time), then the lower triangular matrix was 
matched to the upper triangular matrix. The remaining 19 en- 
semble members are then chosen using the same method, but 
using the control state z (j (i, 1) for i = 1, 2, . . . , / as the central 
reference point. Thus, we have 



z f Q (i,m) = z J Q (i, 1) + 



■ft 



for i = 1,2, 



, / and m = 2, 3, 



(4) 



, 20, where z Q is a matrix 



with ensemble members as columns. Finally, the initial fast 
variables for the forecast are set equal to the initial fast variables 
of the truth, z f Q (j, m) - zg(j, 1) for = /+ 1 , 1+2, and 
m - 1,2, ... , 20. These fast variables are subsequently ignored 
when measuring the shadowing time. 



3.2. Characterizing the ensemble 

Once the ensemble of initial states has been created for each 
hypersphere, the trajectory of each ensemble member is fore- 
cast using the model. Let A be the I x(m- 1) matrix of ensem- 
ble members with the control state at the origin, where m is the 
number of ensemble members (including the control forecast) 
and we are only considering slow variables. If the matrix A acts 
on every vector in the (m - 1 )-dimensional unit sphere, the result 
is a filled-in ellipsoid, provided m > I. The points within the 
ellipsoid can be thought of as representing all potential ensem- 
ble members based upon the known m members (all of which 
lie within this ellipsoid). Although the system is chaotic, over a 
short time span one can approximate the nonlinear evolution of 
the initial sphere of trajectories with this ellipsoid [33]. 

Every two time steps, the ellipsoid is analyzed to deter- 
mine expanding and contracting directions (i.e. those direc- 
tions with positive/negative local finite-time Lyapunov expo- 
nents). This frequency was chosen rather than every time step 
to ensure the identified directions continue to expand/contract. 
Let A = USV T be a singular value decomposition (SVD) of 
A, and let s(f) = S\,S2,-- • ,Si be the singular values (the di- 
agonal entries of S) and U(t) = [uj,U2, ■ • ■ , u„] the matrix of 
left-singular vectors at a given time t. Thus the lengths and di- 
rections of the semi-axes of the ellipsoid approximately encom- 
passing all potential ensemble members are given by s,-U; 1132 1 . 
At each time interval, the vector s(f) is compared to s(t - 1) 




Figure 6: A two-dimensional representation of a 20-member ensemble (/ = 
2). The black point is the control forecast, and the white points are the other 
ensemble members. Let A be the 2x19 matrix of ensemble members with the 
control forecast defined as the origin. The semi-axes are calculated using SVD, 
where s\ and S2 are the two singular values of the matrix A. Note that the SVD 
spectrum can be quite steep for forecasts longer than a few time steps. 



(the singular values from the previous time interval) to deter- 
mine which singular values are decreasing. We then computed 
the dot product between all possible combinations of former 
iii(t - 1) and current u,{t) left singular vectors, with the largest 
dot product chosen to identify direction pairs. Note that on oc- 
casion, a perfect matching was not achieved, but the best guess 
was used. A two-dimensional analog is shown in Fig. [6] Note 
all ensemble members lie within the ellipse, and the control 
state is located at the center by definition of A. SVD allows us 
to identify the ellipse semi-axes s\VL\ and S2VL2- 

An example one-dimensional 40 day ensemble forecast is 
shown in Fig. 0(7 = 6,7 = 16). Note that for the first 20 
days, almost all ensemble members remain close to the truth. 
However, by day 40 the spread of ensemble members covers the 
entire state space, and the ensemble mean no longer accurately 
represents the truth. In addition, Fig.|7]illustrates that the model 
achieves greater extreme values than the system. With less en- 
ergy being dissipated through the fast variables, more energy 
remains oscillating between the slow variables in the model, as 
compared to the system. Note that the reduced regularity of the 
model relative to the system can also be observed. 



4. Forecasting 

Forecasting often fails when the state space dynamics of the 
physical system expand along a previously contracting dimen- 
sion. A two-dimensional illustration of this concept is shown in 
Fig. [8] Because of the uncertainty associated with the initial lo- 
cation of the truth, an /-dimensional sphere of radius cr is used 
to represent the truth. It is from within this sphere that the initial 
ensemble members are chosen. The ellipse representing the en- 
semble members diverges from the trajectory of the truth when 
the system quickly turns in a new direction. In the illustration, 
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Figure 9: Schematic of forecasting with inflation in two dimensions. The black 
line represents the trajectory of the truth. Z ( ' is the ellipsoid encompassing the 
ensemble members at time t, and the dashed line is the trajectory of the ensem- 
ble mean. At each time step, the ellipsoid is inflated by ip along contracting 
directions, where in this case there is only one such direction. Note that fore- 
casting now achieves its goal in contrast to Fig. [8] 



Figure 7: A 20-member ensemble forecast for xj initialized using Eq. {4). The 
black line represents the truth, an integration of the system. The gray lines are 
individual ensemble members, integrations of the model (/ = 6, / = 16). 




Figure 8: Schematic of forecasting in two dimensions. The black line rep- 
resents the trajectory of the truth (z D ). Z, is the ellipsoid encompassing the 
ensemble members at time t, and the dashed line is the trajectory of the ensem- 
ble mean. For nonlinear systems, this ellipsoid is only a good approximation 
of the collections of forecasts for small t. Note that the time steps 0, 1 , 2, 3 are 
an incomplete sampling of all time steps. The goal of a forecast is to achieve 
an ensemble which has a nonempty intersection with i a F . Here forecasting fails 
after time 2. 



forecasting has failed at time F as there is no overlap between 
the /-dimensional sphere around the truth z a F and the ensem- 
ble ellipsoid Z F . For all experiments, the number of expanding 
directions observed throughout the forecast was non-constant. 
This unstable dimension variability is a cause of shadowing 
failures well-documented in the literature 

4.1. Forecasting with inflation 

As described in the introduction, inflation is performed in 
an attempt to improve the ensemble forecast by expanding the 
ellipsoid along contracting dimensions of the state space. For 
each direction u, in which the ensemble is contracting, ensem- 
ble members are inflated along the / th semi-minor axis of the 
ellipsoid s,u,. Although this adds artificial uncertainty to the 
forecast, Danforth and Yorke argue that the amount is minimal 
1 2211 . If the ellipsoid continues to contract in these dimensions, 
the synthetic uncertainty continuously decreases. However, if 



the system begins to expand along a previously contracting di- 
mension, the inflated ensemble will capture some of the change. 
This concept is illustrated for two dimensions in Fig. [9] 

As described above, the directions of the semi-axes are cal- 
culated using SVD throughout the forecast at periodic time in- 
tervals and s(f) is compared to s(f - 1). For each i in which 
Si < Sj(t - 1), all ensemble members are inflated in direction 
u,-. It is important to note here that if inflation were performed 
in expanding dimensions as well, forecasts would degrade more 
rapidly. By constraining inflation to the contracting dimensions, 
we are minimizing the damage caused in hyperbolic regions of 
state space. 

Let <p be the inflation amount, expressed as a percentage of 
the semi-axis length, and let A be the matrix of ensemble mem- 
bers with the control state at the origin, prior to inflation. Define 



M = I + <p(UiuJ) 



(5) 



where I is the /-dimensional identity matrix, uj is the transpose 
of u,, and let A' = MA. Adding the control state to each vector 
in A' then yields ensemble members inflated by ip in direction 
u,. 

The benefit of inflation for a particular forecast is illustrated 
in Fig.[10]for a two-dimensional projection of an ensemble fore- 
cast (/ = 6). As the trajectory of z" begins to turn clockwise, 
the forecasts (gray lines) diverge from the black line. As repre- 
sented by the gray lines, forecasting fails at location 1 without 
inflation (Fig.fTOh). when there is no longer overlap between the 
cr-sphere surrounding z a and the ensemble ellipsoid. However, 
in Fig.fTOb, the gray lines are inflated at location 1 in previously 
contracting directions (as inferred from SVD of the model en- 
semble). Note that the gaps in the figure represent both an infla- 
tion and a time step. As the figure illustrates, these shifts help 
capture the change in direction for z" . It is important to note 
that no measurements of z fl are used to determine the directions 
in which to inflate. A second inflation for the gray lines oc- 
curs four time steps later at location 2 in panel b. At first there 
is more uncertainty in the ellipsoid encompassing the ensem- 
ble members. However, as z a continues in the same direction, 
the uncertainty is dampened by the contracting dynamics. Note 
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let z" represent only the slow variables for simplicity. The root 
mean square error (RMSE) and anomaly correlation (AC) were 
calculated at each time step as follows: 



(a) No inflation 




(b) Using inflation 

Figure 10: Successful inflation in an ensemble forecast. The black line is the 
true trajectory z", and the gray lines are ensemble member forecasts without in- 
flation (panel a) and with inflation having the potential to occur at a minimum 
of every two time steps (panel b). All trajectories are moving in a clockwise di- 
rection. At location 1, forecasting fails without inflation. Location 2 represents 
a second inflation of the ensemble. Inflation successfully modifies the ensemble 
to include the truth as a plausible member. Note that no observations of z" are 
used during the inflation at locations 1 and 2 (/ = 6, J = 16). 



that although the ellipsoid is analyzed every two time steps, the 
model did not inflate between locations 1 and 2. Also, inflation 
acts symmetrically on an ensemble, leaving the ensemble mean 
unchanged during a single step. However, once the ensemble is 
integrated forward in time, the inflated and non-inflated ensem- 
ble mean forecasts diverge. 

4.2. Results 

For each of the 500 initial hyperspheres, a 50-day model fore- 
cast for the 20 ensemble members was created with and with- 
out inflation. At each point in time, let z* be the mean of the 
20 ensemble members' slow variables. For each hypersphere, 
z* was compared to the slow variables of the truth, where we 



RMSE 



z "lb 



AC 



(z* - z) ■ (z° - z) 
||z*-z|| 2 ||z«-z|| 2 



(6) 



(7) 



where z represents the vector of system climatological averages 
for each slow variable. RMSE and AC calculations were then 
averaged over all hyperspheres for the duration of each forecast. 
For the system, RMSE as 9 represents saturation, at which point 
there is no similarity between the forecast and the truth. 

For the present discussion, we consider the period for which a 
forecast is useful, namely the time for which AC remains above 
0.6 J7|]. For each state space size, the average useful time was 
calculated using the non-inflated forecasts. This yielded du- 
rations of 6.7,8.7, and 17.0 days, corresponding to / = 4,5, 
and 6, respectively. If the inflated and non-inflated forecasts 
crossed the 0.6 threshold during the same time step, the useful 
times were considered indistinguishable. Inflation was deemed 
to have "succeeded" if the duration for an acceptable fore- 
cast improved by more than 5% of the average useful time for 
<p = forecasts. Similarly, the technique was considered to 
have "failed" if the duration worsened by the same threshold. 
If there was any improvement (greater than one time step) with 
inflation, forecasts were denoted "inflation helped." Any mea- 
surable degradation with inflation was labeled "inflation hurt." 
Results for / = 4,5, and 6, with tp = 0.5%, 1%,2%, and 5%, 
are recorded in Fig. [TT]with totals out of 500. Bars above zero 
represent improvements, while bars below zero represent degra- 
dations. Black bars indicate improvements and degradations of 
more than 5% of the average useful time (i.e. the number of 
successes and failures). Other metrics such as RMSE and the 
trajectory of the best ensemble member were considered, but 
they showed the same general trends. For all experiments de- 
scribed by Fig. [TT] J — 16 as illustrated in Fig. [3] 

Although the sample size is small, Fig.QT]suggests that infla- 
tion has a greater effect for larger 7, with the best results exhib- 
ited for 7 = 6 and <p = 1 %. In fact, the number of hyperspheres 
in which inflation helped was an order of magnitude better for 
7 = 5 than for 7 = 4. The vast majority of inflated 7 = 4 fore- 
casts were identical to non-inflated (<p = 0) forecasts. These 
numbers roughly doubled for 7 = 6 relative to 7 = 5, likely due 
to the increased regularity in z" for 7 = 6 observed in Fig. [3] 
Unfortunately, the same trend is seen for forecasts in which in- 
flation failed. By capturing one of the sudden changes in the 
trajectory of the truth using inflation, the model can better fore- 
cast the truth for a longer duration. Alternatively, this increase 
in the occurrence of inflation helping might be evidence that er- 
ror introduced by inflation in contracting directions has a min- 
imal effect, as argued by Danforth and Yorke I122ll . This effect 
would likely be strongest in systems with a steep singular value 
spectrum, where the expanding dimensions play a more domi- 
nant role. 

Naturally, the success of inflation is also strongly dependent 
on (p, with increasing tp generally corresponding to a larger 
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Figure 11: Inflation forecasting results for 500 hyperspheres. Bars above zero 
show the number of hyperspheres for which inflation improved the forecast, 
while bars below zero show the number of hyperspheres for which inflation 
degraded the forecasts. Black bars indicate improvements and degradations of 
the average shadowing time by more than 5% ("inflation succeeded" and "infla- 
tion failed," respectively). Gray bars represent the categories "inflation helped" 
and "inflation hurt." The number of times inflation improved the forecast gener- 
ally exceeded the number of degradations, both of which tend to increase with 
dimensionality and ip. 



number of both successes and failures. By increasing the in- 
flation amount, the ellipsoid encompassing the ensemble mem- 
bers is more likely to overlap with the trajectory of the truth, and 
thus capture unstable dimension variability events. On the other 
hand, inflation can make a forecast worse. By inflating in di- 
rections incommensurate with the truth, the method introduces 
additional error, which can be quite significant when amplified 
by nonlinearity. The number of hyperspheres in which inflation 
failed also naturally increases with ip. However, as mentioned 
above, provided the contracting directions continue to contract, 
the error introduced is hopefully minimal. 

For a given set of experimental parameters, the averaged use- 
ful time over all 500 hyperspheres is used for analysis. How- 
ever, the ip — forecast provides an individual useful time for 
each hypersphere. Whether or not inflation helps or hurts (i.e. 
reduces or increases the individual useful time) is not dependent 
upon the magnitude of the individual useful time. For example, 
inflation is equally likely to help (or hurt) an individual forecast 
whose non-inflated useful time was 3 days or 20 days. 

As shown in Fig. Q~T] for many of the hyperspheres infla- 
tion successfully improved the duration for which a forecast is 
useful. The lack of harm caused by this technique in a global 
sense can be seen by averaging RMSE and AC over the hyper- 
spheres in which inflation was successful as shown in Fig. [12] 
for 7 = 6,/= 16. Both plots show small improvement rel- 
ative to forecasts made with no inflation. In section VI, we 
demonstrate the ability to isolate conditions likely to result in 
successful inflation as in Fig.flOl 

Conducting experiments for larger dimensional systems was 
difficult due to computation time limitations for estimating the 
local dynamics of the system attractor (see Section III). To re- 




10 15 20 

Time (days) 
(a) RMSE 



35 



U 
< 




10 15 20 25 30 35 
Time (days) 

(b) AC 

Figure 12: Averaged RMSE (a) and AC (b) for the / = 6, J = 16 system. 
The dashed line represents no inflation, while the solid line is with inflation 
of 1%. Figures show the averaged data for all 500 hyperspheres, while insets 
only include the hyperspheres for which inflation succeeded. Given the small 
(positive) change in average forecast statistics, the method is unlikely to harm 
forecasts, and may be adjusted to isolate particular events for which inflation 
can be quite successful (see Section VI). 



duce the computational burden, initial ensemble members were 
chosen using a spherically-symmetric Monte Carlo approach 
(i.e. nearby initial conditions but with no attractor-imposed 
structure) for these larger systems {I = 7, 8,..., 40). The 
I - 6,1 - 16, and 7 = 21 systems had substantially greater av- 
erage useful times than the other systems (by a factor of five). 
These specific systems all exhibit increased regularity in their 
time series (see Fig. [3}, presumably a result of some resonance 
between the slow and fast variables. For these systems the 
model is most sensitive to inflation, with small <p greatly im- 
proving forecasts and larger <p greatly degrading forecasts. The 
less-regular systems had smaller incidence of success and fail- 
ure, with results similar to those observed for the 1-4 and 
7 = 5 systems. However, the use of a spherically-symmetric 
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Figure 13: Schematic of the redefining of a shadowing ensemble in two dimen- 
sions. The star represents the known location of the truth. The solid circles are 
supplemented by a collection of points meant to approximate the overlap of the 
ellipse and the sphere around the truth. The open circles are previous ensemble 
members that will now be redefined within the overlap, which is approximated 
by a new ellipsoid. Note that the x\ , X% axes are exaggerated for visual ease. 

Monte Carlo approach to initialize an ensemble degraded the 
performance of inflation when compared to the more sophisti- 
cated ensemble initialization technique. 

5. Shadowing and Stalking 

5.7. Methods 

Shadowing with an ensemble presents an alternative method 
for assessing the quality of a model's predictions. Contrasting 
the forecasting experiments described previously, in this con- 
text the trajectory of the truth is known a priori to within a given 
uncertainty (cr), and its location at a given time is represented 
by an /-dimensional sphere of radius cr. At regular intervals, the 
intersection between the ellipsoid encompassing the ensemble 
members and the cr-sphere is approximated. The ensemble el- 
lipsoid is then redefined to represent the intersection, while the 
other ensemble members are redefined to fall within the over- 
lap. Provided this intersection is nonempty, the trajectories of 
some ensemble members yield accurate representations of re- 
ality. A schematic of this technique is shown in Fig. [13] Note 
that although the actual truth lies outside the ensemble ellip- 
soid, some ensemble members are within <x. The trajectories 
of these ensemble members can be forecast forward in time, 
and the process is repeated. If some ensemble members remain 
in the intersection for the entire forecast, then their trajectories 
have successfully cr-shadowed the truth. When inflation is ap- 
plied to the ensemble ellipsoid (as described previously), this 
aggressive form of shadowing is called stalking. 

5.2. Results 

The same initial conditions utilized in the forecasting exper- 
iments were used to create 50-day stalking ensembles. How- 
ever, every four time steps the ensemble ellipsoid was redefined 
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Figure 14: Stalking results for 500 hyperspheres. Bar categories are identical 
to those in Fig.[TT] Inflation improves and worsens the ability to shadow with 
similar probability. Counts tend to slightly increase with the dimension of the 
system. 

as described in Fig. Q~3] This frequency was chosen to be long 
enough to ensure that appropriate contracting dimensions could 
be identified. During the definition of the overlap between the 
ensemble ellipsoid and the truth sphere, semi-axes can be al- 
tered drastically, requiring several time steps to reflect the local 
dynamics. The cr-shadowing distance (uncertainty around the 
truth) was chosen to be 10% of the climatological range in each 
dimension, namely twice the spread in the initial ensemble. A 
new collection of 20 ensemble members were then chosen ly- 
ing within the redefined ellipsoid. As before, RMSE and AC 
were calculated and compared for both shadowing and stalk- 
ing experiments. Counts for the number of hyperspheres for 
which inflation succeeded, failed, helped, and hurt as described 
in Section IV were tallied. Results for / = 4,5, and 6, with 
ip = 0.5%, 1%, 2%, and 5%, are given in Fig. FBI 

The shadowing experiments exhibited a greater benefit from 
inflation when the dimensionality / is large, although the differ- 
ence is less than before. However, the number of hyperspheres 
for which inflation succeeded does not seem to trend mono- 
tonically with <p over the range 0.5% - 5%. Since the number 
of hyperspheres in which inflation succeeded roughly equals 
the number for which inflation helped, we can conclude that 
when stalking is constructive, it improves the shadowing time 
by more than one day. On the other hand, the number of hy- 
perspheres in which inflation succeeded is nearly equivalent to 
the number for which inflation failed for all experiments. Thus, 
inflation can also have a negative effect upon shadowing. 

As with forecasting, stalking experiments were extended 
to higher dimensional systems using a spherically-symmetric 
Monte Carlo approach for ensemble initialization. Similar to 
the lower dimensional systems in Fig. [14] counts for the number 
of hyperspheres for which inflation succeeded and failed were 
relatively close. For small values of inflation, there were often 
errors in the procedure by which corresponding ensemble semi- 
axes directions are calculated between time steps, resulting in 
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an increase in the number of failures. With larger inflation, the 
algorithm can better determine corresponding directions during 
redefinition of the ensemble. The highly regular I - 6,1 - 16, 
and 7 = 21 systems were less distinguished for stalking experi- 
ments. Since the truth is known throughout, ensemble members 
are forced by the redefinition to remain nearby regardless of the 
regularity of the attractor. 

6. Targeted Inflation 

Despite taking advantage of local stability estimates, infla- 
tion as defined previously was naive in the sense that it occurred 
in all contracting directions, including those which push fore- 
casts away from the truth. Thus, this naive inflation had the 
potential to both improve and degrade a forecast. The final ob- 
jective of this research was to isolate the instances in which in- 
flation is beneficial to the forecast, and identify the systematic 
tendencies of these situations. With this knowledge, at each 
time step the decision of whether inflation should occur, and in 
which contracting directions, could be based upon the charac- 
teristics of the local state space of the system. For these targeted 
inflation experiments, the 1-5 system was chosen for further 
study; it does not have the same regularity as the 1-6 system 
(Fig. |3j, but is more regular than some of the higher dimen- 
sional systems. 

For a given time step in the forecast, the shape of the attrac- 
tor near the current state z/ was approximated as follows. First, 
the state 50 time steps prior in the forecast z{_ 50 was considered. 
By adding uniform random noise in every direction on the order 
of <x, 1000 analogs located near (within cr) z^_ 5 „ were chosen. 
These 1000 analogs were forecast forward 50 time steps, and 
the covariance matrix associated with z{ was calculated. The 
eigenvectors of this forecast state covariance matrix yield an 
approximation of the shape of the attractor near z/ . Finally, the 
dot product of these eigenvectors and the proposed inflation di- 
rection (determined by SVD to be contracting) was calculated. 
The larger this projection, the more likely the proposed inflation 
direction is consistent with the local shape of the attractor. If the 
resulting dot product was above a set threshold fi and the direc- 
tion was determined to be contracting, then inflation proceeded 
as before. Otherwise, no inflation occurred. For all values of <p, 
fi was set at 0,0.8, and 0.9, reducing the number of inflations 
by roughly 0%, 60%, and 80%, respectively. 

The goal here is to avoid inflating in contracting directions 
likely to harm the forecast. By initializing the analogs ran- 
domly, we are able to sample a new collection of state space 
directions, distinct from those spanned by the inflated ensem- 
ble. 

6.1. Forecasting 

Forecasting results are given in Fig. Q3] Adjusting the num- 
ber of inflations throughout the forecast by changing fi clearly 
improves the method. Note first that results from a cutoff of 
l-i = agree with Fig. QT| The small differences (most no- 
table for <p = 1%) can be attributed to starting from different 
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Figure 15: Targeted inflation for forecasting experiments on 500 hyperspheres. 
Bar categories are the same as those defined in Fig. 1111 Note that /j = 0.9 
substantially reduces the instances of failure, while not significantly reducing 
the instances of success (/ = 5). 

initial points in the attractor. For ip = 2% and tp = 5%, increas- 
ing fi has the effect of decreasing the counts; with fewer infla- 
tions, ensembles more closely resemble the non-inflated fore- 
cast {ip = 0). It is important to note that although the number 
of hyperspheres for which inflation succeeds remains relatively 
constant, the number for which inflation fails is substantially re- 
duced as ft is increased. This trend indicates that the technique 
is successfully determining when inflation will indeed be use- 
ful, and is inflating less frequently when it is likely to worsen a 
forecast. For example, consider tp = 5% in Fig. [15] With naive 
inflation {fi = 0), the technique succeeded for 39 hyperspheres 
and failed 42 times. Setting fi = 0.9 reduced the number of 
successes by 23% to 30, however, it eliminated 7 1 % of the fail- 
ures, leaving only 12. When fi is increased beyond 0.9 though, 
inflation occurs too rarely and the number of hyperspheres for 
which inflation succeeds is greatly reduced (not shown). 

For (p = 0.5% and <p = 1%, the trend is a little more com- 
plicated. Still, fi = 0.9 resulted in a remarkably better ratio 
of number of hyperspheres succeeding versus failing. With 
if = 0.5%, fi = 0.8 surprisingly increased the number of hyper- 
spheres succeeding and helping relative to fi — 0. For tp = 1%, 
fi = 0.8 had notably few hyperspheres for which inflation suc- 
ceeded. This indicates that although experiments were run for 
500 different hyperspheres, there is still variability between tri- 
als based upon initial conditions and chaotic system dynamics. 

6.2. Stalking 

As shown in Fig. [16] the results for stalking experiments are 
more variable. For some algorithm parameters, such as fi — 0.8 
for tp = 0.5% and tp = 1%, targeted inflation appears to suc- 
cessfully increase the number of hyperspheres for which in- 
flation succeeded, while decreasing the number for which in- 
flation failed. However, the same fi for greater inflation had 
smaller differences, even increasing the number of trajectories 
for which inflation failed. One possible explanation for the 
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Figure 16: Targeted inflation for stalking experiments on 500 hyperspheres. 
Bar categories are the same as those defined in Fig. 1111 Results varied, with 
more work needed to optimize the success of inflation for stalking experiments 
(/ = 5). 

mixed results in stalking experiments is the experimental de- 
sign. Redefining the ensemble every fourth time step reduces 
the ability to identify expanding/contracting dimensions. This 
can lead to difficulties in matching ellipsoid semi-axis direc- 
tions from one time step to the next. Adjusting yU did alter the 
results, however more work is needed to optimize the success 
of inflation for stalking experiments. 

7. Discussion 

Using the Lorenz '96 coupled system as an analog for atmo- 
spheric dynamics, we examined the potential for using inflation 
to improve model forecasts of the Earth's atmosphere. In the 
first experiments, forecasts were made with and without infla- 
tion, and the resulting trajectories were compared to the system 
truth. In the second part of this work, the existence of a shad- 
owing trajectory with and without inflation was assessed. The 
trajectory of the truth was known a priori, and only the clos- 
est ensemble members were considered every few time steps. 
In the former, the technique is referred to as inflation, and in 
the latter, stalking. Under idealized conditions, inflation was 
shown to be beneficial for both techniques. The final part of the 
work attempted to isolate the ensemble directions most likely 
to benefit from increased uncertainty using targeted inflation. 

For the Lorenz '96 model, by adjusting the number of slow 
and fast variables, and the degree of coupling between them, 
one could vary between regular and highly chaotic systems; the 
greatest regularity was observed for 1-6 and 7=16. However, 
by increasing the number of fast variables (/), greater regularity 
can be achieved for a fixed number of slow variables. In con- 
trast to hyperbolic systems, the present system exhibits unstable 
dimension variability and fluctuating Lyapunov exponents as is 
common in higher dimensional systems [35]. 

For forecasting experiments, inflation showed the potential 
to increase the length of time for which an ensemble forecast 



represents the verifying truth as a plausible member, as illus- 
trated in Fig. [10] Inflation in directions away from the truth, 
however, can degrade a forecast. As a result, greater inflation 
amounts always correspond to an increase in the number of hy- 
perspheres for which inflation failed. For regular systems such 
as I = 6 and I = 16, smaller inflation amounts can often suc- 
cessfully improve a forecast, while the error introduced from 
larger inflation amounts is magnified. However, for less regu- 
lar systems, larger inflation was necessary to improve forecasts. 
For these systems, the uncertainty introduced was less likely to 
have a harmful effect. 

Stalking experiments were used to establish the existence of 
a model trajectory that remains close to the truth for a given 
time, and assess how inflation affects this trajectory. For these 
experiments, results were relatively constant for systems with 
greater than six dimensions. Inflation was shown to be success- 
ful at about the same rate at which it hindered shadowing. Thus, 
the utility of stalking is dependent upon the ability to isolate the 
hyperspheres for which inflation helped. 

Finally, by quantifying the local dynamics at each time step 
and restricting inflation to directions consistent with the lo- 
cal attractor, the instances for which inflation helped were 
isolated with moderate success for forecasting experiments. 
The number of hyperspheres for which inflation failed could 
be decreased without significantly decreasing the number for 
which inflation succeeded. Stalking experiments with targeted 
inflation yielded mixed results, with additional modifications 
needed to attain the same degree of success. In practice how- 
ever, the truth is not known a priori. Thus, the improvements 
observed with targeted inflation for forecasting experiments are 
encouraging for application in NWP. 

When shadowing physical systems such as the Earth's at- 
mosphere, one is presented with the challenges of uncertainty 
in the initial condition, sensitive dependence, and model er- 
ror. Model error is estimated to dominate the forecast error 
in weather systems for the first three days ||r7ll . Nevertheless, 
modelers must attempt to provide the best predictions possible 
given the circumstances. This research demonstrates that in- 
flation has the potential to aid in our attempts to model chaotic 
physical systems, provided we are able to isolate the local shape 
of the attractor and tune the inflation parameters successfully. 
Targeted inflation is most likely to be beneficial for highly non- 
linear events (e.g. hurricanes), where sudden changes in the 
state space trajectory are prevalent. Further exploration of the 
design parameters of targeted inflation are likely to result in in- 
creased effectiveness. 

One potential application of this technique to NWP is in 
maintaining the spread of a forecast ensemble. For unbiased 
ensembles, an effective spread is equal to the RMSE of the en- 
semble mean with respect to the analysis [36], but NWP en- 
sembles typically fail to retain this ideal beyond the first few 
days. Of the various means by which one could increase the 
ensemble spread artificially, targeted inflation preserves the en- 
semble orientation relative to the local model attractor and has 
the potential to improve the ensemble mean. 

The major disadvantage of targeted inflation, relative to the 
naive version, is the frequently required computation of the 
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shape of the attractor. For NWP models, where the number of 
degrees of freedom is (9(10 10 ), this computation would be pro- 
hibitive. However, the operational forecasting centers compute 
new analyses and ensemble forecasts roughly every six hours. 
Thus, the required information for an estimate of the attractor 
shape would be available cost free in terms of computational 
time. 
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